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Abstract 

We present an exploratory study of a one-parameter family of covariant, 
non-perturbative lattice gauge-fixing conditions, that can be implemented 
through a simple Monte Carlo algorithm. We demonstrate that at the 
numerical level the procedure is feasible, and as a first application we 
examine the gauge dependence of the gluon propagator. 
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1 Introduction 



In recent years, hadron spectroscopy has become the most popular field of ap- 
plication of lattice QCD. Hadron masses, decay constants and matrix elements 
for semileptonic decays are routinely computed from Green functions of com- 
posite hadron fields and expectation values of current operators between hadron 
states. However, the predictive power of lattice QCD is not limited to these 
kinds of calculations. In particular. Green functions of individual quark and 
gluon fields can be computed as well as expectation values of operators between 
quark and gluon states. The interest of this approach is twofold. Firstly, one can 
non-perturbatively compute renormalisation constants for composite operators 
by sandwiching them between quark states 0. These calculations form a crucial 
step in extracting physical information from lattice calculations. Secondly, quark 
and gluon Green functions are interesting objects in their own right. Being the 
most fundamental computable quantities in QCD, they are expected to contain 
direct information on the mechanism of colour confinement and chiral symmetry 
breaking. Also, they allow a direct determination from first principles of the 
running QCD coupling and may be relevant for understanding the physics of 
pomeron exchange from the point of view of QCD [Q. 

Unlike hadronic operators. Green functions for quark and gluons must be 
defined in a fixed gauge. For most of the applications described above, it is a 
crucial problem to disentangle gauge dependent features from gauge independent 
ones. For example, if a dynamical mass were extracted from a nonperturbative 
study of the gluon propagator, a physical interpretation may be attached to it 
only if one can obtain reasonable evidence that such a mass does not depend on 
the gauge chosen, at least within a class of gauges. For this reason, it would be of 
great interest to be able to define and implement on the lattice a whole family of 
nonperturbative gauge conditions, by varying continously some gauge parameter. 

In this paper we describe a numerical study of one such class of gauges. In 
section 2 we recall the formulation of the gauge condition, in the framework of 
the Feynman path integral, and the corresponding Monte Carlo algorithm. In 
Section 3 we discuss the numerical performance of the algorithm and we present 
a preliminary study of the gauge dependence of the gluon propagator. Finally, 
In Section 4 we sketch our agenda for the future. 

2 Global Gauge-Fixing 
2.1 General Framework 

We start from Wilson's lattice gauge model, defined by the gauge invariant par- 
tition function: 




(1) 
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where /3 = for a gauge group SU{N) and Sw is the standard Wilson action. 
The formula for the expectation value of an observable O is: 

< O >w= J dU e-^^^[^] 0[U]. (2) 

It is well known (Elitzur's theorem) that if 0[U] is a local, gauge dependent 
function the above expression vanishes. 

In [0 a general procedure for nonperturbative lattice gauge-fixing was pro- 
posed. They defined a modified partition function by simply inserting a factor 1 
in (0): 

Zmod ^JdU e-'^^-I^l r'[U] I dg e-P'^'^^^'\ (3) 

where the dg integration runs over the group of lattice gauge transformations. 
F\p^] is a generic function of the links such that it is not invariant under general 
gauge transformations and I\U] is defined as 

I[U] = jdg e-/5^^'^[^^l. (4) 

Clearly I\U] is gauge invariant. As usual, the gauge-transformed link is defined 
as 

Ul{n)^g{n)U,{n)g\n + ^^). (5) 

Formula (|^) corresponds to the first step of the standard Faddeev-Popov 
gauge-fixing procedure. However, unlike what would happen in the continuum, 
Zmod = Zw is a finite quantity because the group of gauge transformations on a 
finite lattice is compact. Thus Zmod can provide a new definition for the expec- 
tation value of 0[U]: 

< O >_,^ Z-l, J dU e-^^-[^l r'[U] J dg e"^^^^^!^^! 0[U% (6) 

In general, if 0[U] is gauge dependent its expectation value in the modified scheme 
does not vanish. If 0[U] is gauge invariant then < O >w=< O >mod, so that (|^) 
defines a consistent, nonperturbative gauge-fixing procedure. 
By defining 

<0[U]>G=I-'[U] j dg e-^^'"^^'''^ 0{U% (7) 
< O >mod can be cast in the form: 

/rff/e-^^^[^] <0\U]>G ^,rr, 

< O >mod= ^ jdUe-^sMu\ =< < 0[U] >c >w • (8) 

The above expression indicates that in the gauge-fixed model the expectation 
value of a gauge dependent quantity 0[U] is obtained in two steps. First one 
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associates with 0\U] the gauge invariant function < 0[U] >g, which has the 
form of a Gibbs average of 0[f/^] over the group of gauge transformations, with 
a statistical weight factor e~^^^ fiu^]^ Then one takes the average of < 0[U] >g 
a la Wilson. 

This suggests the following numerical algorithm: 

1. generate a set of link configurations Ui,...Un, weighted by the Wilson 
action, via the usual gauge invariant Monte Carlo algorithm for some value 
of/9; 

2. use each of the Ui as a set of quenched "bonds" in a new Monte Carlo 
process, where the dynamical variables are the local gauge group elements 
g{n) located on the lattice sites. These are coupled through the hnks Ui, 
according to the effective Hamiltonian F[Uf]. In this way one can produce 
for every link configuration Ui an ensemble of gauge-related configurations, 
weighted by the Boltzmann factor exp(-/3M^F[f/f ]). We call < 0[Ui] >g 
the average of a gauge dependent observable O with respect to such an 
ensemble, in the spirit of Formula (^. 

3. Finally, the expectation value < O >mod is simply obtained from the Wilson 
average of the < 0[t/i] >g, i-e.: 

1 ^ 

< O >mod- T7 E < Om >G . (9) 

i=l 

In the above scheme can be interpreted as a gauge parameter, which 
determines the effective temperature 1//3M^ of the Monte Carlo on the group 
of gauge transformations. In the following we will refer to such a Monte Carlo 
process as GFMC. The goal of our project is to simulate the system for many 
different values of M^, corresponding to different gauge choices. This would 
enable us to study the gauge dependence of relevant lattice quantities. 

2.2 A Convenient Class of Gauges 

For our numerical study we will use the gauge-fixing function 

F[U^]^-Y: ReTr{Ul{n)), (10) 

where the sum runs over all lattice links. There are many motivations for such a 
choice. Recalling (|^), it turns out that as an effective Hamiltonian for the GFMC 
process the above function describes a simple nearest-neighbour interaction of 
the variables g{n), with couplings given by the corresponding (quenched) links. 
The g{n) can be interpreted as S't/(3)-valued spins, and the couplings are also 
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elements of SU{3). In this sense ( ]T0| ) defines a classical, 4- dimensional SU{3) 
spin-glass with SU (3) couplings. 

For small values of (3 and M^, numerical results may be compared with ana- 
lytical ones, derived from a strong coupling expansion (see for details). Also, 
a weak coupling expansion of a continuum version of this model has been per- 
formed by Fachin 0]. The latter, which is valid for any value of M^, may prove 
extremely useful when investigating the continuum limit of our lattice system. 

Besides practical advantages, the choice (^) also has a theoretical interest 
as the continuum version of this class of gauges was proposed as a possible 
solution to the Gribov problem in the Landau gauge . Landau gauge-fixing has 
similar features on the lattice and in the continuum formulation of the theory. In 
particular, it has been shown that Gribov copies also exist on the lattice Q]. The 
connection of such a scheme with the Landau gauge appears when one studies 
the behaviour of (0) as a function of the g variables, for fixed U. Then it turns 
out that the stationary points of F[U^] correspond to link configurations that 
satisfy the lattice version of the Landau gauge condition. All such configurations 
correspond to Gribov copies. In particular, those corresponding to local minima 
of F[U^] also satisfy a positivity condition for the lattice Faddeev-Popov operator 
0. As a consequence, in the limit — oo, the above gauge- fixing is equivalent 
to the so-called minimal Landau gauge condition, which prescribes to pick up on 
every gauge orbit the field configuration corresponding to the absolute minimum 
of F[U^] |T^. We will not discuss the Gribov problem in further detail, as it is 
not central to our present purpose. Here our main goal is to determine whether 
the above scheme for lattice gauge-fixing can be efficiently simulated numerically. 



3 Numerical Results 

3.1 Performance of the Algorithm and Thermodynamics 

For our exploratory study we considered quenched QCD on a 8^ lattice at /3 = 5.7. 
All the numerical work was performed on single-processor Alpha workstations, 
located at the University of Edinburgh. We first generated link configurations, 
weighted by the Wilson action, using a hybrid-overrelaxed algorithm, where both 
the Cabibbo-Marinari (CM) pseudo-heat-bath and overrelaxed (OR) updates 
were performed on three SU{2) subgroups. Next, as described in the previous 
section, for each link configuration we produced an ensemble of gauge-related con- 
figurations, weighted by the Boltzmann factor exp{—l3 M'^ F[U^]) . This was done 
for many different values of jSM'^ and again the GFMC sweep was a combination 
of CM and OR updates. 

Before analysing the above system we studied a simpler one, for testing pur- 
poses. This was obtained by setting all the link variables to the identity, cor- 
responding to the limit P ^ oo, and then generating "pure gauge" configura- 
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tions, weighted by exp{— M'^F[P]). The resulting system has the form of a 
four-dimensional SU{3) spin model with ferromagnetic couplings. We examined 
its thermodynamics in view of a comparison with the three-dimensional case, 



which has been studied in the literature llll. We found evidence for a first order 



phase transition, at = 1/T = 0.635. This was obtained by studying the 
M^— dependence of the specific heat of the system, defined as 

where E = F[P] and the gauge group average is defined according to (0). In the 
three-dimensional first order transition was also observed by Kogut and 

collaborators. 

Turning to the system at /5 = 5.7, we analysed a range of values for 
up to jSM"^ = 2.4. The first problem that we addressed was the study of the 
thermalisation of the GFMC process. In fact, by recalling the analogy of F[U^] 
with a spin-glass Hamiltonian, one expects metastable states to appear at low 
temperature (i.e. for large M^). From the point of view of the algorithm one has 
to make sure that, for a chosen range of values of M^, the stochastic process can 
efficiently visit all the states and does not get trapped in a metastable one. To this 
end we chose one of the link configurations and performed several GFMC runs on 
it, using different seeds for the random number generator. For each value of 
at least four different seeds where used. We then analysed the evolution in GFMC 
time of two observables. One was the expectation value of F[U^], corresponding 
to the average energy of the system < E[U] >g (cfr. Eq. (^), the other being 
the zero-momentum gluon 2-point function, whose precise definition will be given 
in the following subsection. The important point here is that while the former 
quantity is a local one, the latter depends on the dynamics of long-wavelength 
modes of the system. For this reason we used the 2-point function to measure 
the autocorrelation time of our algorithm. When plotting the evolution of the 
observables against the number of sweeps, metastable states could be identified 
as the occurrence of "false" plateaux, i.e. long-lived stable values that eventually 
"decay" into the real vacuum. One example is shown in Fig]^, which is a plot 
of the zero momentum gluon 2-point function vs. the number of sweeps. At this 
stage, by changing the number of CM and OR updates in the GFMC sweep, we 
were able to optimise the performance of the algorithm at every value of M^, 
in order to get rid of metastable states. In general we found that by increasing 
the ratio of OR vs. CM moves the algorithm became more efficient. In the high 
temperature region, /3M^ < 0.8, a GFMC sweep composed of one CM sweep and 
two OR updates resulted in an autocorrelation time of the order of one GFMC 
sweep. For /?M^ > 0.8 we had to increase the number of OR updates, finally 
setting the GFMC sweep to be a combination of 1 CM and 10 OR updates. This 
resulted in an autocorrelation time of the order of 4 sweeps. 



5 



20 



% 

? ^ A <^ ^^^^ . , <<to A. A^- 



15 r ^<?.^ 



10 — 



5 — 











3000 



1000 2000 
sweeps 

Figure 1: zero momentum gluon 2-point function vs. GFMC sweeps in a typical 
thermalisation run. 



In summary, we were able to tune the algorithm so as to obtain reasonable 
evidence that for pM^ < 2.4 the GFMC process could thermalise correctly. A 
similar analysis for other link configurations confirmed the pattern found with the 
one used for tuning. We then performed a detailed study of the thermodynamics. 
As expected from the behaviour of the algorithm, two separate regions in the 
parameter space could be identified, corresponding to a strong coupling regime 
for /?M^ < 0.8, and a weak coupling one, for /?M^ > 0.8. A phase transition 
seems to separate them. The nature of the transition was analysed again in terms 
of the specific heat, defined now as 

C^od = jj^ d(iTM^ " ^ " ^ ^ ^^^^ 

(see FigJ^). Notice that in (|12D the average is also taken with respect to link 
configurations (cfr. Eq. (^). The transition still appeared to be first order, as 
in the ferromagnetic case, and seemed to occur for the same value of the critical 
temperature. 

In the weak coupling region our numerical data are obtained from 21 link 
configurations. For each of them an ensemble of 60 gauge-related configurations 
was generated. In the strong coupling region the statistics is smaller. Statistical 
errors, obtained from a jackknife analysis, are shown in Figs. 2-5. However, in 
most cases such error is negligible. 

Fig.H shows our numerical results for < E >mod vs. /3M^ and the correspond- 
ing analytical result from a strong coupling expansion up to next to leading order 
[If. The agreement is perfect up to /?M^ ^ 0.7. 
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Figure 2: specific heat vs. j3 M^. 

3.2 Gauge Dependence of the Gluon Propagator 

Having gained confidence that we could correctly generate ensembles of gauge- 
related configurations for a significant range of values, we attempted a study 
of the gluon propagator. This quantity has been extensively studied at the non- 
perturbative level both in the continuum, mainly through Schwinger-Dyson equa- 
tions |jl2|, and numerically on the lattice [|l^ . 



One issue of particular relevance is the behaviour of the propagator for mo- 
menta p {). It has been advocated that nonperturbative effects dynamically 
generate a gluon mass which removes the infrared pole of the perturbative prop- 
agator. More generally, a behaviour softer than a pole has been advocated by 
many authors. Lattice studies seem to provide some support for the mass gener- 
ation hypothesis |T^, but one important point is the possible gauge dependence 



of such a mass. Within the limits of a preliminary investigation, we attempted to 
gain some insight into this issue by computing the propagator for several values 
of M2. 

We recall the lattice definition of the gluon field in terms of the link variables 
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M^) = " - ^Tr ( ) , (13) 



2iagQ 3 \^ 2iagQ 

where a is the lattice spacing and go is the bare coupling constant. By Fourier 
transforming the above field, one can define the bare lattice n-point gluon Green 
functions, in momentum space: 

<^S2.../.„(^'l'^'2, ...,Pn) = {Af,,{pi)Af,,{p2) . . . Af,^{pn))mod, (14) 
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Figure 3: < E >mod vs. P M^. The dotted line is the result from the strong 
coupling expansion. 

where {■)mod obviously indicates the average according to (^) and momentum 
conservation implies pi + p2 + • • • + Pn = 0. 
The gluon propagator is then defined as 

G^;^lM = {A,,{p)A,,i-p))^^od. (15) 
In particular, we studied the following quantity: 

G^Lrip) = E {A,ip)A,i-p))^,,. (16) 

As discussed in the previous section, a comparison of our results for finite 
values of with those obtained in the minimal Landau gauge was of particular 
interest, as this gauge corresponds to the —>■ oo limit of our scheme. Strictly 
speaking, a numerical implementation of the minimal Landau gauge is not feasible 
on the lattice as it would require finding the absolute minimum of a spin-glass-like 
Hamiltonian. In practice, we assumed that the role of multiple minima could be 
neglected, i.e. we identified the minimal Landau gauge with the gauge obtained 
by imposing that F[U^] attained a local minimum. Such an approximation has 



been widely used in the literature [13 



In Fig.|^ we plot GfJ^i^^{p) vs. p in lattice units, for a range of values of /3M^ 
and in the minimal Landau gauge. Our data indicate a strong dependence on the 
gauge parameter M^. In particular, for the lowest value of /3M^, which sits on 
the phase transition, the propagator is essentially momentum independent, being 
dominated by a large effective mass generated by thermal fluctuations. 
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Figure 4: gluon propagator vs. momentum, for various values of /9M^, and in 
the minimal Landau gauge. 



Focusing on the behaviour at p = 0, we show in Fig.|^ GsJ^iariP = 0) as a 
function of /?M^. The dotted line indicates the Landau gauge value. Again, 
the M^— dependence appears quite substantial and strongly correlated with the 
phase transition, as expected. G^^J^i^j.{p = 0) appears to increase monotonically 
with towards the minimal Landau gauge asymptotic value. 

Before drawing any firm conclusion from our results we need to address the 
issue of the existence of the continuum limit. We plan to investigate this point in 
the near future by repeating our calculation for different values of (3. In particular, 
by defining a renormalised gluon propagator as in , one should be able to check 
that such a quantity does not depend on /?, for fixed M^. 

In this connection, it should be observed that at the numerical level one 
does not expect our scheme to survive the continuum limit for all values of M^. 
On one hand, since in the limit — > the gauge-fixing effect disappears, 
we expect the signal-to-noise ratio for a gauge dependent quantity to get worse 
with decreasing M^, with the average value eventually approaching zero. This is 
consistent with the data in Figs. 3-5. Obviously one does not expect the scheme 
to have a continuum limit for values of jSM"^ in the strong coupling region, where 
in fact the gluon propagator is suppressed by thermal fluctuations (see Figs. 4-5). 
On the other hand, for large values of we may face increasing difficulties in 
the thermalisation process, so that in practice one hopes that a significant range 
of values of will remain accessible when increasing (3. 
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Figure 5: gluon propagator at p = vs. /?M^. 

4 Conclusions 

We have performed a preliminary numerical study of a nonperturbative lattice 
gauge-fixing scheme, where a whole class of gauges can be implemented by vary- 
ing a gauge parameter. An efficient algorithm was devised and tested against 
analytical results in the strong coupling region. As a first application we studied 
the gauge dependence of the gluon propagator, which appears to be quite sub- 
stantial, in particular at zero momentum. This preliminary result, if confirmed 
in a more complete analysis, would suggest that no physical meaning could be 
attached to a dynamically generated gluon mass. In order to establish the rele- 
vance of our results to continuum physics, we plan to increase the range of values 
for the gauge-fixing parameter and to repeat our calculation for higher values of 
(3 and larger lattices. Such a study will be presented in a future publication. 
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